This vignette illustrates the use of INLA for spatial prediction using examples from Blangiardo and Cameletti (2015) and Illian, Sørbye, and Rue (2012). For prediction of continuous spatial processes, the Lindgren, Rue, and Lindström (2011) stochastic partial differential equations (SPDE) approach is used to approximate the process through an areal Gaussian Markov random field (GMRF) representation. Finally, Log-Gaussian Cox process models are fit using the pseudodata approach of Simpson et al. (2016).
Blangiardo and Cameletti (2015) section 6.1.
Toy dataset from Blangiardo and Cameletti (2015). A spatial process is observed via a random sampling plan that concentrates observations in the lower left of the unit square.
# Plot the data.
plot(s2 ~ s1, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0), data = SPDEtoy, pch = 19, asp = 1, main = 'Toy Data')
# Create a mesh for the SPDE method and then plot it.
toy_mesh <- inla.mesh.2d(as.matrix(SPDEtoy[,c('s1', 's2')]), max.edge = c(0.1, 0.2))
plot(toy_mesh, asp = 1)
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
# SPDE projector matrix for estimation.
A_est <- inla.spde.make.A(toy_mesh, as.matrix(SPDEtoy[,c('s1', 's2')]))
# Initialize exponential covariance structure for SPDE.
spde <- inla.spde2.matern(mesh = toy_mesh, alpha = 2)
# Set up stack for estimation.
stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = spde$n.spde)
stack_est <- inla.stack(data = list(y = SPDEtoy$y), A = list(A_est), effects = list(c(stack_index, list(intercept = 1))), tag = 'est')
# Create a grid for prediction.
toy_nx <- 50
toy_ny <- 50
toy_grid <- expand.grid(x = seq(0, 1, length.out = toy_nx), y = seq(0, 1, length.out = toy_ny))
# SPDE projector matrix for prediction.
A_pred <- inla.spde.make.A(mesh = toy_mesh, loc = as.matrix(toy_grid))
# Set up stacks for prediction.
stack_latent <- inla.stack(data = list(xi = NA), A = list(A_pred), effects = list(stack_index), tag = 'pred_latent')
stack_response <- inla.stack(data = list(y = NA), A = list(A_pred), effects = list(c(stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
stacks <- inla.stack(stack_est, stack_latent, stack_response)
# Fit the model with INLA.
toy_fit <- inla(
y ~ -1 + intercept + f(spatial_field, model = spde),
data = inla.stack.data(stacks),
control.predictor = list(A = inla.stack.A(stacks), compute = TRUE)
)
# Output posterior summaries.
toy_fit$summary.fixed
toy_fit$summary.hyperpar
# Extract posterior mean of latent spatial field.
index_latent <- inla.stack.index(stacks, tag = 'pred_latent')$data
post_mean <- toy_fit$summary.linear.predictor[index_latent, 'mean']
post_sd <- toy_fit$summary.linear.predictor[index_latent, 'sd']
# Plot the posterior mean and SD of the latent spatial field.
plot(im(matrix(post_mean, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior Mean of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
plot(im(matrix(post_sd, nrow = toy_nx, ncol = toy_ny), xrange = range(toy_grid$x), yrange = range(toy_grid$y)), main = 'Posterior SD of Spatial Field')
points(SPDEtoy$s1, SPDEtoy$s2, col = rgb(SPDEtoy$y / max(SPDEtoy$y), 0, 0, 0.5), pch = 20)
Example from Møller and Waagepetersen (2007), Beilschmiedia pendula Lauraceae locations in a plot in Panama. bei dataset in spatstat (Baddeley and Turner 2005). The illustrations below fit a stationary LGCP model with no covariates, a Matèrn covariance function with \(\nu = 1\), and INLA’s default (flat?) priors for all parameters.
# Plot the full point pattern.
plot(bei, pch = '.', cols = 'black', main = 'Realized Point Pattern')
# Take a sample of quadrats and plot the observed point pattern.
set.seed(84323)
N_QUADS <- 10
QUAD_SIZE <- 50
w_edge <- Frame(bei)$xrange[1]
e_edge <- Frame(bei)$xrange[2]
s_edge <- Frame(bei)$yrange[1]
n_edge <- Frame(bei)$yrange[2]
botleft <- cbind(
runif(N_QUADS, w_edge, e_edge - QUAD_SIZE),
runif(N_QUADS, s_edge, n_edge - QUAD_SIZE)
)
bei_interior <- lapply(seq_len(nrow(botleft)), function(r){return(
cbind(
botleft[r, 1] + c(0, 0, QUAD_SIZE, QUAD_SIZE),
botleft[r, 2] + c(0, QUAD_SIZE, QUAD_SIZE, 0)
)
)})
bei_win <- do.call(
union.owin,
apply(botleft, 1, function(x){return(
owin(x[1] + c(0, QUAD_SIZE), x[2] + c(0, QUAD_SIZE))
)})
)
bei_hole <- bei[complement.owin(bei_win, frame = Frame(bei))]
bei_samp <- bei[bei_win]
bei_window_full <- Window(bei)
plot(bei_hole, main = 'Observed Region with Holes', pch = '.', cols = 'black')
plot(bei_window_full, main = 'Observed Sample')
plot(bei_win, add = TRUE)
plot(bei_samp, pch = '.', cols = 'black', add = TRUE)
The mesh should be fine in observed regions to accurately represent complex local structure but can be coarsened in unobserved regions where the model cannot infer as much detail. I am including margins to explore model behavior away from the observed regions.
# Parameters to experiment with.
MAX_EDGE_LENGTH <- 25
MAX_EDGE_EXT <- 50
MARGIN <- 100
# Mesh covering the site.
bei_boundary <- inla.mesh.segment(loc = do.call(cbind, vertices.owin(Window(bei))))
bei_full_mesh <- inla.mesh.create(
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_full_spde <- inla.spde2.matern(mesh = bei_full_mesh)
plot(bei_full_mesh, asp = 1, main = 'Fine Mesh')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh')
# Mesh including a margin outside the site.
margin_mesh <- inla.mesh.2d(
loc = bei_full_mesh$loc[,1:2], # Include nodes from site.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_spde <- inla.spde2.matern(mesh = margin_mesh)
plot(margin_mesh, asp = 1, main = 'Fine Mesh with Coarse Margin')
points(bei, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Margin')
# Meshs with coarser resolution in quadrats.
quad_hole <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = bei_interior[[i]], grp = i - 1))
})
)
bei_hole_mesh0 <- inla.mesh.create(
boundary = list(bei_boundary, quad_hole),
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_hole0_spde <- inla.spde2.matern(mesh = bei_hole_mesh0)
plot(bei_hole_mesh0, asp = 1, main = 'Fine Mesh with Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Holes')
bei_hole_mesh <- inla.mesh.create(
loc = bei_hole_mesh0$loc[,1:2], # Include nodes from mesh with holes.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_hole_spde <- inla.spde2.matern(mesh = bei_hole_mesh)
plot(bei_hole_mesh, asp = 1, main = 'Fine Mesh with Coarse Holes')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh with Coarse Holes')
# Meshs with finer resolution in quadrats.
quad_bnd <- do.call(
inla.mesh.segment,
lapply(seq_along(bei_interior), function(i){
return(inla.mesh.segment(loc = apply(bei_interior[[i]], 2, rev), grp = i - 1))
})
)
bei_samp_mesh0 <- inla.mesh.create(
boundary = quad_bnd,
refine = list(max.edge = MAX_EDGE_LENGTH)
)
bei_samp0_spde <- inla.spde2.matern(mesh = bei_samp_mesh0)
plot(bei_samp_mesh0, asp = 1, main = 'Fine Mesh in Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Mesh in Quadrats')
bei_samp_mesh <- inla.mesh.create(
loc = bei_samp_mesh0$loc[,1:2], # Include nodes from mesh in quads.
boundary = bei_boundary,
refine = list(max.edge = MAX_EDGE_EXT) # Fill in the rest with a coarser triangulation.
)
bei_samp_spde <- inla.spde2.matern(mesh = bei_samp_mesh)
plot(bei_samp_mesh, asp = 1, main = 'Coarse Mesh with Fine Quadrats')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Coarse Mesh with Fine Quadrats')
# Meshes with varying resolution in quadrats and a margin.
margin_hole <- inla.mesh.2d(
loc = bei_hole_mesh$loc[,1:2], # Include nodes from mesh with holes.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_hole_spde <- inla.spde2.matern(mesh = margin_hole)
plot(margin_hole, asp = 1, main = 'Coarse Holes and Margin')
points(bei_hole, pch = 19, cex = 0.25, col = 'red')
title('Coarse Holes and Margin')
margin_samp <- inla.mesh.2d(
loc = bei_samp_mesh$loc[,1:2], # Include nodes from quads.
offset = MARGIN,
max.edge = MAX_EDGE_EXT # Fill in the rest with a coarser triangulation.
)
margin_samp_spde <- inla.spde2.matern(mesh = margin_samp)
plot(margin_samp, asp = 1, main = 'Fine Quadrats and Coarse Margin')
points(bei_samp, pch = 19, cex = 0.25, col = 'red')
title('Fine Quadrats and Coarse Margin')
# Identify which mesh nodes are in the oberserved region
# and create SPDE projectors for spatial mapping.
NPIX_X <- 400
NPIX_Y <- 200
obs_full <- rep(0, margin_mesh$n)
obs_full[inla.over_sp_mesh(as(Window(bei), 'SpatialPolygons'), margin_mesh, 'vertex')] <- 1
proj_margin_mesh <- inla.mesh.projector(margin_mesh, dims = c(NPIX_X, NPIX_Y))
obs_hole <- rep(0, margin_hole$n)
obs_hole[inla.over_sp_mesh(as(Window(bei_hole), 'SpatialPolygons'), margin_hole, 'vertex')] <- 1
proj_margin_hole <- inla.mesh.projector(margin_hole, dims = c(NPIX_X, NPIX_Y))
obs_samp <- rep(0, margin_samp$n)
obs_samp[inla.over_sp_mesh(as(Window(bei_samp), 'SpatialPolygons'), margin_samp, 'vertex')] <- 1
proj_margin_samp <- inla.mesh.projector(margin_samp, dims = c(NPIX_X, NPIX_Y))
Count the points in grid cells and fit a Poisson GLMM with the observed area in the cell as an exposure variable. The plots of the cell coutns are blank (white) where cells have no observed area.
NGRID_X <- 40
NGRID_Y <- 20
centers <- gridcenters(
dilation(bei_window_full, max(NGRID_X, NGRID_Y)),
NGRID_X, NGRID_Y)
dx <- sum(unique(centers$x)[1:2] * c(-1, 1)) / 2
dy <- sum(unique(centers$y)[1:2] * c(-1, 1)) / 2
bei_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_df))){
bei_df$count[r] <- sum(bei$x >= bei_df$x[r] - dx &
bei$x < bei_df$x[r] + dx &
bei$y >= bei_df$y[r] - dy &
bei$y < bei_df$y[r] + dy)
bei_df$area[r] <- area(Window(bei)[owin(c(bei_df$x[r] - dx, bei_df$x[r] + dx), c(bei_df$y[r] - dy, bei_df$y[r] + dy))])
}
)
user system elapsed
0.456 0.004 0.463
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_df$area > 0, bei_df$count, NA), nrow = length(unique(bei_df$x)))), unique(bei_df$x), unique(bei_df$y), unitname = 'meters'), ncolcours = range(bei_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(bei_window_full, border = 'white', add = TRUE)
points(bei, pch = '.', col = 'black')
# SPDE projector matrix for estimation.
full_A_est <- inla.spde.make.A(
margin_mesh,
as.matrix(bei_df[bei_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
full_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_spde$n.spde)
full_stack_est <- inla.stack(data = list(count = bei_df$count[bei_df$area > 0], larea = log(bei_df$area[bei_df$area > 0])), A = list(full_A_est), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
full_A_pred <- inla.spde.make.A(mesh = margin_mesh, loc = as.matrix(bei_df[,c('x', 'y')]))
# Set up stacks for prediction.
full_stack_latent <- inla.stack(data = list(xi = NA), A = list(full_A_pred), effects = list(full_stack_index), tag = 'pred_latent')
full_stack_response <- inla.stack(data = list(count = NA), A = list(full_A_pred), effects = list(c(full_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
full_stacks <- inla.stack(full_stack_est, full_stack_latent, full_stack_response)
# Fit the model with INLA.
system.time(
bei_full_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(full_stacks),
control.predictor = list(A = inla.stack.A(full_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
896.228 0.404 121.059
# Output posterior summaries.
bei_full_fit$summary.fixed
bei_full_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, bei_full_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, bei_full_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
bei_hole_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_hole_df))){
bei_hole_df$count[r] <- sum(bei_hole$x >= bei_hole_df$x[r] - dx &
bei_hole$x < bei_hole_df$x[r] + dx &
bei_hole$y >= bei_hole_df$y[r] - dy &
bei_hole$y < bei_hole_df$y[r] + dy)
bei_hole_df$area[r] <- area(Window(bei_hole)[owin(c(bei_hole_df$x[r] - dx, bei_hole_df$x[r] + dx), c(bei_hole_df$y[r] - dy, bei_hole_df$y[r] + dy))])
}
)
user system elapsed
1.148 0.004 1.150
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_hole_df$area > 0, bei_hole_df$count, NA), nrow = length(unique(bei_hole_df$x)))), unique(bei_hole_df$x), unique(bei_hole_df$y), unitname = 'meters'), ncolcours = range(bei_hole_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
hole_A_est <- inla.spde.make.A(
margin_hole,
as.matrix(bei_hole_df[bei_hole_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
hole_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_hole_spde$n.spde)
hole_stack_est <- inla.stack(data = list(count = bei_hole_df$count[bei_hole_df$area > 0], larea = log(bei_hole_df$area[bei_hole_df$area > 0])), A = list(hole_A_est), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
hole_A_pred <- inla.spde.make.A(mesh = margin_hole, loc = as.matrix(bei_hole_df[,c('x', 'y')]))
# Set up stacks for prediction.
hole_stack_latent <- inla.stack(data = list(xi = NA), A = list(hole_A_pred), effects = list(hole_stack_index), tag = 'pred_latent')
hole_stack_response <- inla.stack(data = list(count = NA), A = list(hole_A_pred), effects = list(c(hole_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
hole_stacks <- inla.stack(hole_stack_est, hole_stack_latent, hole_stack_response)
# Fit the model with INLA.
system.time(
bei_hole_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_hole_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(hole_stacks),
control.predictor = list(A = inla.stack.A(hole_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
848.172 0.448 115.200
# Output posterior summaries.
bei_hole_fit$summary.fixed
bei_hole_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, bei_hole_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, bei_hole_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
bei_samp_df <- data.frame(x = centers$x, y = centers$y,
count = NA_integer_, area = NA_real_)
system.time(
for(r in seq_len(nrow(bei_samp_df))){
bei_samp_df$count[r] <- sum(bei_samp$x >= bei_samp_df$x[r] - dx &
bei_samp$x < bei_samp_df$x[r] + dx &
bei_samp$y >= bei_samp_df$y[r] - dy &
bei_samp$y < bei_samp_df$y[r] + dy)
bei_samp_df$area[r] <- area(Window(bei_samp)[owin(c(bei_samp_df$x[r] - dx, bei_samp_df$x[r] + dx), c(bei_samp_df$y[r] - dy, bei_samp_df$y[r] + dy))])
}
)
user system elapsed
0.916 0.004 0.919
par(mar = c(0.5, 0, 2, 2))
plot(im(t(matrix(ifelse(bei_samp_df$area > 0, bei_samp_df$count, NA), nrow = length(unique(bei_samp_df$x)))), unique(bei_samp_df$x), unique(bei_samp_df$y), unitname = 'meters'), ncolcours = range(bei_samp_df$count) %*% c(-1, 1) + 1, main = 'Binned Tree Counts')
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = '#00000040')
# SPDE projector matrix for estimation.
samp_A_est <- inla.spde.make.A(
margin_samp,
as.matrix(bei_samp_df[bei_samp_df$area > 0, c('x', 'y')])
)
# Set up stack for estimation.
samp_stack_index <- inla.spde.make.index(name = 'spatial_field', n.spde = margin_samp_spde$n.spde)
samp_stack_est <- inla.stack(data = list(count = bei_samp_df$count[bei_samp_df$area > 0], larea = log(bei_samp_df$area[bei_samp_df$area > 0])), A = list(samp_A_est), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'est')
# SPDE projector matrix for prediction.
samp_A_pred <- inla.spde.make.A(mesh = margin_samp, loc = as.matrix(bei_samp_df[,c('x', 'y')]))
# Set up stacks for prediction.
samp_stack_latent <- inla.stack(data = list(xi = NA), A = list(samp_A_pred), effects = list(samp_stack_index), tag = 'pred_latent')
samp_stack_response <- inla.stack(data = list(count = NA), A = list(samp_A_pred), effects = list(c(samp_stack_index, list(intercept = 1))), tag = 'pred_response')
# Join all three stacks.
samp_stacks <- inla.stack(samp_stack_est, samp_stack_latent, samp_stack_response)
# Fit the model with INLA.
system.time(
bei_samp_fit <- inla(
count ~ -1 + intercept + f(spatial_field, model = margin_samp_spde),
offset = larea, family = 'poisson',
data = inla.stack.data(samp_stacks),
control.predictor = list(A = inla.stack.A(samp_stacks), compute = TRUE),
verbose = TRUE
)
)
user system elapsed
263.856 0.252 36.574
# Output posterior summaries.
bei_samp_fit$summary.fixed
bei_samp_fit$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, bei_samp_fit$summary.random$spatial_field$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, bei_samp_fit$summary.random$spatial_field$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
This method relies upon the Lindgren, Rue, and Lindström (2011) approximation of the latent Gaussian field as a linear combination of a finite number of basis functions represented as a GMRF on the nodes of a triangulation of the space. Simpson et al. (2016) use the triangulation for numerical integration of the intensity function and show that the LGCP likelihood factors into the joint distribution of independent Poisson random variables corresponding to the points of the point pattern and the nodes of the triangulation. The model fitting proceeds using INLA to fit a Poisson model to pseudodata.
The pseudodata are constructed as follows.
Then \(y_{i} \sim Poisson(\alpha_{i}\eta_{i})\) where \(\log(\eta_{i})\) is the SPDE representation of the GF at the location of the \(i\)th pseudodatum. See the paper for tedious notation regarding the definition of \(\eta_{i}\). Ultimately, the nodes become Poisson random variables with means equal to the intensity at that their respective locations, observed points become Poisson random variables with means of 1, and the likelihood is approximately proportional to
\[\prod_{i=1}^{p+n} \eta_{i}^{y_{i}} \exp(-\alpha_{i} \eta_{i}).\]
(Is there a missing \(\alpha_{i}\)?)
Sampling effort is accomodated by scaling the integration weights by the probability that a point would have been observed at that node (given the sampling plan), i.e. nodes in observed regions have the \(\alpha_{i}\) defined above and nodes in unobserved regions have \(\alpha_{i} = 0\). Weights can be scaled by values other than 0 or 1 to account for e.g. distance sampling with a (known) detection function (Yuan et al. 2017) or a (known) false negative rate for detection equipment. This adjustment assumes the point process of interest is observed through a known thinning process and then allows inference back to the intensity function of the unthinned process.
full_pts <- cbind(bei$x, bei$y)
# Contruct the SPDE A matrix for nodes and points.
full_nV <- margin_mesh$n
full_nData <- dim(full_pts)[1]
full_LocationMatrix <- inla.mesh.project(margin_mesh, full_pts)$A
full_IntegrationMatrix <- sparseMatrix(i = 1:full_nV, j = 1:full_nV, x = rep(1, full_nV))
full_ObservationMatrix <- rbind(full_IntegrationMatrix, full_LocationMatrix)
# Get the integration weights.
full_IntegrationWeights <- diag(inla.mesh.fem(margin_mesh)$c0)
full_E_point_process <- c(obs_full * full_IntegrationWeights, rep(0, full_nData))
# Create the psuedodata.
full_fake_data <- c(rep(0, full_nV), rep(1, full_nData))
# Fit model to full site.
full_formula <- y ~ -1 + intercept + f(idx, model = margin_spde) # No covariates.
full_data <- list(y = full_fake_data, idx = 1:full_nV, intercept = rep(1, full_nV))
system.time(
result_full <- inla(
formula = full_formula,
data = full_data,
family = 'poisson',
control.predictor = list(A = full_ObservationMatrix),
E = full_E_point_process,
verbose = TRUE
)
)
user system elapsed
416.864 0.420 60.643
result_full$summary.fixed
result_full$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_mesh, result_full$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
hole_pts <- cbind(bei_hole$x, bei_hole$y)
# Contruct the SPDE A matrix for nodes and points.
hole_nV <- margin_hole$n
hole_nData <- dim(hole_pts)[1]
hole_LocationMatrix <- inla.mesh.project(margin_hole, hole_pts)$A
hole_IntegrationMatrix <- sparseMatrix(i = 1:hole_nV, j = 1:hole_nV, x = rep(1, hole_nV))
hole_ObservationMatrix <- rbind(hole_IntegrationMatrix, hole_LocationMatrix)
# Get the integration weights.
hole_IntegrationWeights <- diag(inla.mesh.fem(margin_hole)$c0)
hole_E_point_process <- c(obs_hole * hole_IntegrationWeights, rep(0, hole_nData))
# Create the psuedodata.
hole_fake_data <- c(rep(0, hole_nV), rep(1, hole_nData))
# Fit model to site with holes.
hole_formula <- y ~ -1 + intercept + f(idx, model = margin_hole_spde) # No covariates.
hole_data <- list(y = hole_fake_data, idx = 1:hole_nV, intercept = rep(1, hole_nV))
system.time(
result_hole <- inla(
formula = hole_formula,
data = hole_data,
family = 'poisson',
control.predictor = list(A = hole_ObservationMatrix),
E = hole_E_point_process,
verbose = TRUE
)
)
user system elapsed
356.252 0.452 51.817
result_hole$summary.fixed
result_hole$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_hole, result_hole$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
samp_pts <- cbind(bei_samp$x, bei_samp$y)
# Contruct the SPDE A matrix for nodes and points.
samp_nV <- margin_samp$n
samp_nData <- dim(samp_pts)[1]
samp_LocationMatrix <- inla.mesh.project(margin_samp, samp_pts)$A
samp_IntegrationMatrix <- sparseMatrix(i = 1:samp_nV, j = 1:samp_nV, x = rep(1, samp_nV))
samp_ObservationMatrix <- rbind(samp_IntegrationMatrix, samp_LocationMatrix)
# Get the integration weights.
samp_IntegrationWeights <- diag(inla.mesh.fem(margin_samp)$c0)
samp_E_point_process <- c(obs_samp * samp_IntegrationWeights, rep(0, samp_nData))
# Create the psuedodata.
samp_fake_data <- c(rep(0, samp_nV), rep(1, samp_nData))
# Fit model to quadrat-sampled site.
samp_formula <- y ~ -1 + intercept + f(idx, model = margin_samp_spde) # No covariates.
samp_data <- list(y = samp_fake_data, idx = 1:samp_nV, intercept = rep(1, samp_nV))
system.time(
result_samp <- inla(
formula = samp_formula,
data = samp_data,
family = 'poisson',
control.predictor = list(A = samp_ObservationMatrix),
E = samp_E_point_process,
verbose = TRUE
)
)
user system elapsed
293.136 0.676 44.242
result_samp$summary.fixed
result_samp$summary.hyperpar
# Plot surface.
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx$mean)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
par(mar = c(0.5, 0, 2, 2))
plot(im(t(inla.mesh.project(proj_margin_samp, result_samp$summary.random$idx$sd)),
xrange = Frame(bei)$x + c(-MARGIN, MARGIN),
yrange = Frame(bei)$y + c(-MARGIN, MARGIN),
unitname = c('meter', 'meters')),
main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
inlabruinlabru has a wrapper function to fit an LGCP with INLA, with a relatively easy-to-use interface for defining models and predicting arbitrary functions of latent variables. However, it is poorly documented, slow, and the documentation seems to imply that it does not support variable sampling effort (even though this appears to work).
bei_full_spdf <- as.SpatialPoints.ppp(bei)
cmp_full <- coordinates ~ mySmooth(map = coordinates, model = margin_spde) + Intercept
system.time(
bei_full_lgcp <- lgcp(cmp_full, bei_full_spdf, E = obs_full, options = list(verbose = TRUE))
)
user system elapsed
3417.044 1.320 445.024
bei_full_lgcp$summary.fixed
bei_full_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_full <- predict(bei_full_lgcp, pixels(margin_mesh), ~ mySmooth)
plot(lambda_full, main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
plot(lambda_full['sd'] ,main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
points(bei, pch = '.', col = 'white')
bei_hole_spdf <- as.SpatialPoints.ppp(bei_hole)
cmp_hole <- coordinates ~ mySmooth(map = coordinates, model = margin_hole_spde) + Intercept
system.time(
bei_hole_lgcp <- lgcp(cmp_hole, bei_hole_spdf, E = obs_hole, options = list(verbose = TRUE))
)
user system elapsed
2955.588 1.184 387.283
bei_hole_lgcp$summary.fixed
bei_hole_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_hole <- predict(bei_hole_lgcp, pixels(margin_hole), ~ mySmooth)
plot(lambda_hole, main = 'Posterior Mean of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
plot(lambda_hole['sd'], main = 'Posterior SD of Spatial Field')
plot(Window(bei_hole), border = 'white', add = TRUE)
points(bei_hole, pch = '.', col = 'white')
bei_samp_spdf <- as.SpatialPoints.ppp(bei_samp)
cmp_samp <- coordinates ~ mySmooth(map = coordinates, model = margin_samp_spde) + Intercept
system.time(
bei_samp_lgcp <- lgcp(cmp_samp, bei_samp_spdf, E = obs_samp, options = list(verbose = TRUE))
)
user system elapsed
789.304 1.484 112.311
bei_samp_lgcp$summary.fixed
bei_samp_lgcp$summary.hyperpar
# Plot posterior means and posterior sd.
lambda_samp <- predict(bei_samp_lgcp, pixels(margin_samp), ~ mySmooth)
plot(lambda_samp, main = 'Posterior Mean of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
plot(lambda_samp['sd'], main = 'Posterior SD of Spatial Field')
plot(Window(bei), border = 'white', add = TRUE)
plot(Window(bei_samp), border = 'white', add = TRUE)
points(bei_samp, pch = '.', col = 'white')
All three methods give similar results for the full dataset and the dataset with holes, even when the gridding method uses a coarse grid. The intercept and random effects are shifted from method to method, but theses are not separately indentifyable and the shifts cancel each other out. The methods each have different artifacts and edge effects apparent in the results from the sampled dataset. The pseudodata approach is the fastest except when a very coarse grid is used and a small region is observed.
Baddeley, Adrian, and Rolf Turner. 2005. “Spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.
Blangiardo, Marta, and Michela Cameletti. 2015. Spatial and Spatio-Temporal Bayesian Models with R-INLA. Wiley.
Illian, Janine B, Sigrunn H Sørbye, and Håvard Rue. 2012. “A Toolbox for Fitting Complex Spatial Point Process Models Using Integrated Nested Laplace Approximation (Inla).” The Annals of Applied Statistics, 1499–1530.
Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (4): 423–98.
Møller, J, and RP Waagepetersen. 2007. “Modern Spatial Point Process Modelling and Inference.” Scandinavian Journal of Statistics 34: 643–711.
Simpson, Daniel, Janine B Illian, Finn Lindgren, Sigrunn H Sørbye, and Havard Rue. 2016. “Going Off Grid: Computationally Efficient Inference for Log-Gaussian Cox Processes.” Biometrika 103 (1): 49–70.
Yuan, Yuan, Fabian E Bachl, Finn Lindgren, David L Borchers, Janine B Illian, Stephen T Buckland, Haavard Rue, Tim Gerrodette, and others. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” The Annals of Applied Statistics 11 (4): 2270–97.